Merge Desjardins and Ashton trees

Author

Claudia Zirión-Martínez

Published

February 10, 2025

Setup

Libraries

Code

library(RRphylo)
library(manipulate)
library(ape)
library(phytools)
library(ggtree)
library(tidyverse)
library(RColorBrewer)
library(ggnewscale)
library(patchwork)
source("scripts/metadata_colors.R")

Paths

Code

metadata_ashton_desj_all_weavepop_H99 <- 
    "data/processed/metadata_ashton_desj_all_weavepop_H99.csv"

desj_tree_path <- 
    "data/raw/CryptoDiversity_Desjardins_Tree.tre"
desj_tree_out_path <- 
    "data/processed/tree_desjardins.newick"
desj_tree_out_plot <- 
    "results/trees/tree_desjardins.png"

ashton_tree_path <- 
    "data/raw/2017.06.09.all_ours_and_desj.snp_sites.mod.fa.cln.tree"
ashton_metadata_path <- 
    "../Crypto_Ashton/config/metadata_all_ashton_and_vni_desj.csv"
ashton_tree_out_path <- 
    "data/processed/tree_ashton.newick"
ashton_tree_unrooted_plot <- 
    "results/trees/tree_ashton_unrooted.png"
ashton_tree_rooted_plot_pdf <- 
    "results/trees/tree_ashton.pdf"
ashton_tree_out_plot <- 
    "results/trees/tree_ashton.png"


merged_tree_out_path <- 
    "data/processed/tree_merged.newick"
merged_tree_branchlengths_plot <- 
    "results/trees/tree_merged_branchlengths.png"
merged_tree_plot <- 
    "results/trees/tree_merged.png"
merged_tree_small_plot <- 
    "results/trees/tree_merged_small.png"

Metadata

Use the metadata table that has all the samples included in the final Crypto_Desjardins_Ashton dataset and H99 (n = 1056).

Code

metadata <- read.delim(
    metadata_ashton_desj_all_weavepop_H99,
    header=TRUE,
    sep=",")
summary <- metadata %>% 
    group_by(dataset, lineage) %>% 
    summarize(count = n())
summary

dataset	lineage	count
Ashton	VNI	668
Desjardins	VNBI	122
Desjardins	VNBII	64
Desjardins	VNI	185
Desjardins	VNII	16
Reference	VNI	1

Make separate dataframes for each metadata field.

Code

metadata$vni_subdivision <- factor(metadata$vni_subdivision,
                            levels = names(sublineage_colors))
metadata$country_of_origin <- factor(metadata$country_of_origin,
                                levels = names(country_colors))

sublineage <- metadata %>%
                filter(lineage == "VNI")%>%
                select(strain, vni_subdivision)%>%
                column_to_rownames("strain")%>%
                droplevels()
lineage <- metadata %>%
            select(strain, lineage)%>%
            column_to_rownames("strain")
dataset <- metadata %>%
            select(strain, dataset)%>%
            column_to_rownames("strain")
source <- metadata %>%
            select(strain, source)%>%
            column_to_rownames("strain")
country <- metadata %>%
            select(strain, country_of_origin)%>%
            column_to_rownames("strain")

Desjardins tree

Import the raw Desjardins tree

Code

desj_tree <- read.tree(desj_tree_path)

Reroot the tree at the middle of the branch leading to VNII

Code

VNII_root <- getMRCA(desj_tree, c("C2","C12"))
edge_length <- subset(desj_tree$edge.length, desj_tree$edge[,2] == VNII_root)
desj_tree <- reroot(desj_tree, VNII_root, edge_length/2)

Write rooted Desjardins tree

Code

write.tree(desj_tree, file = desj_tree_out_path)

Keep only the names of the countries in the Desjardins dataset to have a proper legend.

Code

country_desj <- levels(droplevels(country[rownames(country) %in% desj_tree$tip.label, ]))

Code

d <- ggtree(desj_tree, layout = "circular", size = 0.1) +  
            geom_tiplab(aes(label = label),
                        size = 0.8, align =TRUE, 
                        linetype = "dashed", linesize = .03)+
            geom_treescale(x=-0.03, y=0, width=0.01, fontsize =2, offset=30)

d1 <- gheatmap(d, lineage, width=.08, colnames=FALSE, offset=0.05) +
            scale_fill_manual(values = lineage_colors, name="Lineage",
                                na.translate = FALSE)+
            guides(fill = guide_legend(order = 1))+
            new_scale_fill()

d2 <- gheatmap(d1, sublineage, width=.08, colnames=FALSE, offset=0.08) +
            scale_fill_manual(values = sublineage_colors, 
                            name="VNI Sublineage", na.translate = FALSE, 
                            limits = levels(sublineage$vni_subdivision))+
            guides(fill = guide_legend(order = 2))+
            new_scale_fill()

d3 <- gheatmap(d2, source, width=.08, colnames=FALSE, offset=0.11) +
        scale_fill_manual(values = source_colors, name="Source", na.translate = FALSE)+
        guides(fill = guide_legend(order = 3))+
        new_scale_fill()

d4 <- gheatmap(d3, country, width =.08, colnames=FALSE,offset=0.14)+
    scale_fill_manual(values=country_colors, 
                      name="Country", na.translate=TRUE,
                      limits = country_desj)+
    labs(title = "Desjardins tree")+
    guides(fill = guide_legend(order =4))+
        theme(legend.position = "right",
                legend.direction = "vertical",
                legend.title = element_text(size=9),
                legend.text=element_text(size=7),
                legend.key.size = unit(0.3, "cm"),
                plot.margin = margin(0, 0, 0, 0, "cm"))
d4

Ashton tree

Import the raw Ashton tree

Code

ashton_tree_unrooted <- read.tree(ashton_tree_path)

Rename tips of the Desjardins samples because they have SRA run accessions in this tree but strain names in the other one

Code

ashton_tree_unrooted$tip.label <- sapply(ashton_tree_unrooted$tip.label, function(x) {
    if (x %in% metadata$run) {
        metadata$strain[metadata$run == x]
    } else {
        x
    }
})

Get the samples that are present in the tree but absent from the metadata of the final dataset

Code

tips_missing_from_final_dataset <- setdiff(ashton_tree_unrooted$tip.label, metadata$strain)

Compare the list of strains missing from metadata with the oringinal Ashton metadata

Code

ashton_metadata <-read.delim(
    ashton_metadata_path,
    header=TRUE, sep=",")
samples_missing_from_dataset <- ashton_metadata %>%
    filter(strain %in% tips_missing_from_final_dataset)%>%
    select(sample, strain, lineage, VNI_subdivision)
samples_missing_from_dataset

sample	strain	lineage	VNI_subdivision
ERS542414	15277_3#7	VNI	VNIa-4
ERS542415	15277_3#8	VNI	VNIa-4
ERS542595	15277_3#45	VNI	VNIa-4
ERS542403	15277_3#1	VNI	VNIa-4
ERS542456	15277_3#18	VNI	VNIa-4
ERS542410	15277_3#5	VNI	VNIa-5
ERS542411	15277_3#6	VNI	VNIa-5
	CNS_1465	VNI	VNIa-93
ERS542584	15277_3#42	VNI	VNIa-93
ERS542502	14893_1#16	VNI	VNIa-93

The CNS_1465 strain was not available for download and the rest had bad quality alignments.

Root Ashton tree at the middle of the branch leading to VNIa

Code

VNIa_root <- getMRCA(ashton_tree_unrooted, c("AD3-95a","Tu259-1"))
edge_length <- subset(ashton_tree_unrooted$edge.length, 
    ashton_tree_unrooted$edge[,2] == VNIa_root)
ashton_tree <- reroot(ashton_tree_unrooted, VNIa_root, edge_length/2)

Write rooted Ashton tree

Code

write.tree(ashton_tree, file = ashton_tree_out_path)

Code

pu <- ggtree(ashton_tree_unrooted, layout = "circular", size = 0.1) +
            geom_tippoint(aes(color = sublineage[match(label, rownames(sublineage)), "vni_subdivision"]),
                            size = 2)+
            scale_color_manual(values = sublineage_colors, 
                                name="VNI Sublineage", 
                                na.translate = FALSE,
                                limits = levels(sublineage$vni_subdivision))+
            labs(title = "Unrooted")+
            theme(legend.position = "none",
                plot.title = element_text(hjust = 0.5))

p <- ggtree(ashton_tree, layout = "circular", size =0.1)+
            geom_tippoint(aes(color = sublineage[match(label, rownames(sublineage)), "vni_subdivision"]),
                            size = 2)+
            scale_color_manual(values = sublineage_colors, 
                                name="VNI Sublineage", 
                                na.translate = FALSE,
                                limits = levels(sublineage$vni_subdivision))+
        labs(title = "Rooted")+
        theme(legend.position = "right",
                plot.title = element_text(hjust = 0.5),
                legend.direction = "vertical",
                legend.title = element_text(size=11),
                legend.text=element_text(size=9),
                legend.key.size = unit(0.5, "cm"))

pu | p

Code

m <- ggtree(ashton_tree, layout = "circular", size = 0.1) +  
            geom_tiplab(aes(label = label), size = 0.5, align =TRUE, 
                    linetype = "dashed", linesize = .03)+
            geom_treescale(x=-0.03, y=0, width=0.01, fontsize =2, offset=30)

m1 <- gheatmap(m, dataset, width=.05, colnames=FALSE, offset=.025) +
            scale_fill_manual(values = dataset_colors, name="Dataset", na.translate = FALSE)+
            guides(fill = guide_legend(order = 1))+
            new_scale_fill()

m2 <- gheatmap(m1, sublineage, width=.05, colnames=FALSE, offset=.035) +
            scale_fill_manual(values = sublineage_colors, name="VNI Sublineage", na.translate = FALSE,
                                limits = levels(sublineage$vni_subdivision))+
            guides(fill = guide_legend(order = 2))+
            new_scale_fill()

m3 <- gheatmap(m2, source, width=.05, colnames=FALSE, offset=0.045) +
        scale_fill_manual(values = source_colors, name="Source", na.translate = FALSE)+
        guides(fill = guide_legend(order = 3))+
        labs(title = "Ashton tree")+
        new_scale_fill()

m4 <- gheatmap(m3, country, width =.05, colnames=FALSE,offset=0.055)+
    scale_fill_manual(values=country_colors, 
                      name="Country", na.translate=FALSE,
                      limits = levels(country$country_of_origin))+
    guides(fill = guide_legend(order =4))+
        theme(legend.position = "right",
                legend.direction = "vertical",
                legend.title = element_text(size=9),
                legend.text=element_text(size=7),
                legend.key.size = unit(0.3, "cm"),
                plot.margin = margin(0, 0, 0, 0, "cm"))

m4